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ABSTRACT 


T temperature, k 


This paper describes the first phase of an T* 

effort to develop multidimensional models of 
Stirling engine components; the ultimate goal is to t 
model an entire engine working space. More specif- 
ically, this paper describes parallel plate and u 

tubular heat exchanger models with emphasis on the 
central part of the channel (i.e., ignoring hydro- u* 
dynamic and thermal end effects). The model 
assumes: laminar, incompressible flow with con- u av 

stant thermophysical properties. In addition, a 
constant axial temperature gradient is imposed. u 0 

The governing equations, describing the model, have 
been solved using Crack-Nicloson f ini te-di fference x 
scheme. Model predictions have been compared with 
analytical solutions for oscillating/reversing y 

flow and heat transfer in order to check numerical 
accuracy. The simplifying assumptions will later 
be relaxed to permit modeling of compressible, a e 

laminar/turbulent flow that occurs in Stirling heat 
exchanger. 


Excellent agreement has been obtained for the 
model predictions with analytical solutions avail- 
able for both flow in circular tubes and between 
parallel plates. Also the heat transfer computa- 
tional results are in good agreement with the heat 
transfer analytical results for parallel plates. 

NOMENCLATURE 

dp/dx axial pressure gradient, N/m 3 


normalized temperature = T/ (X Ax) 
time, s 

axial velocity, m/s 
normalized axial velocity = u/u 0 
cross-stream-averaged velocity, m/s 
character i sti c axial velocity, m/s 
axial distance, m 

cross stream distance, measured from 
plate centerl ine, m 

effective averaged thermal diffusivity, 
“the thermal diffusivity that will 
produce the same axial heat flux as 
obtained from radial conduction and axial 
convection heat transfer" = 


2 

Tu dy dt, m /s 


normalized effective averaged thermal 
diffusivity = a e /(u> Ax 2 ) 


U) 

2irX(R 0 - R.) 


f time-averaged Fanning friction factor = 



f s unidirectional Fanning friction factor 

Pr Prandtl number 

p pressure, N/m 2 

Re Reynolds number = 2 u av (R 0 - Rp/v 

Ri half width of the channel, m 


Ax tidal displacement "the cross-stream- 

averaged maximum axial distance which the 
fluid elements travel during one half 
period of osci 1 lation" * 


r ir / a* *0 



Ax* normalized tidal displacement * Ax/(u 0 /u) 
9T 

» X axial temperature gradient, k/m 


R 0 half distance between two consecutive 

plates, m 

r radial location, measured from channel 

centerl ine, m 


p fluid dynamic viscosity, kg/m/s 
v fluid kinematic viscosity, m 2 /s 
p density, kg/m 3 


1 





M « 


t w wall shear stress, N/m 2 

4> pressure gradient phase angle, (deg) 

- 0 >t( 360/2-rr) 

w frequency, r ad/s 

w* Valensi number * [w(R 0 - R,) 2 ]/v 

INTRODUCTION 

NASA is developing free-piston Stirling 
engines for space power generation. These engines 
are designed with computer codes which simulate the 
oscillating/reversing flow and oscillating pressure 
level using one-dimensional flow analyses. Confi- 
dence is lacking in the resulting characterizations 
of the various interrelated thermodynamic losses, 
for example, the use of steady-flow friction fac- 
tor and heat transfer correlations throughout the 
engine working space is of questionable validity. 


(1) The flew is laminar, incompressible and 
has constant thermophysical properties. 

iZ) Since we are only concerned with the fluid 
flow in the central part (in the <~direction) of 
the channel, hydrodynamic and thermal end effects 
can be ignored. This assumption can be satisfied 
by using a long channel compared to its width. 

(3) As a consequence to assumption (2) all 
convection terms in the momentum and energy equa- 
tions can be ignored, except for the pu(3T/3x) 
term in the energy equation; since it is the only 
convective driving mechanism for the heat transfer. 

(4) The axial temperature gradient is con- 
stant. This assumption is consistent with the 
experimental finding for a similar setup C2] where 
it was found that the temperature varies linearly 
with x in almost 80 percent of the central length 
of the channel . 


It has not been possible to make accurate 
in-engine measurements of flows, velocities, tem- 
peratures and pressure drops. However, several 
osci 1 lating-f low rigs are being used to learn more 
about engine f 1 ui d mechani cs and heat transfer 
via dynamic measurements. Multidimensional model- 
ing is another process for achieving better under- 
standing of fluid flow and heat transfer phenomena 
in Stirling engine working spaces. Such models can 
be used to study velocity and temperature fields 
within the engine components ; these fields can 
then be compared with the fields implied by use of 
the one-dimensional correlations. Heat transfer 
and fluid friction can also be calculated and com- 
pared with correlation assumptions used in one- 
dimensional codes. It may be possible to validate 
portions of multidimensional models directly with 
data from the oscillating flow rigs. Cleveland 
State University, with support from NASA Lewis 
Research Center, is working on development of 
multidimensional models of Stirling engine compo- 
nents. The first stage of this work, reported 
here, is development of a two-dimensional model of 
parallel plate and tubular heat exchangers. 

One problem encountered by use of any numeri- 
cal code is the evaluation of the errors in the 
solution due to the inherently approximate numeri- 
cal techniques. A practical approach to evaluating 
such errors is to compare the numerical solution 
for a restricted version of the problem to an 
exact analytical solution. This paper reports on 
comparison of parallel plate heat exchanger model 
predictions with an exact analytical solution 
developed by Kurzweg [13. Kurzweg's solution is 
for oscillating flow between parallel plates with 
heat transfer. 


ANALYSIS 


(5) Under the oscillating flow conditions, 
heat transfer is controlled by radial conduction 
and axial convection. Therefore the axial heat 
conduction can be ignored. 

Governing Equations 


Applying the above assumptions to the general 
form of the Navier Stokes equations yields the fol- 
lowing set of equations: 


Continui ty : 


x-Momentum: 


3u 

3x 


3u dp 1 3 / i 

p at " " 37 + j 3?V 


Energy : 

3T l( 3T 13 (S u 3T\ 
p 3t = pU 3x * i 3rV Pr 3rJ 


(1) 


( 2 ) 


(3) 


where i = 1 for circular tubes 

* 0 for flow between parallel plates 


The code developed is capable of predicting both 
circular tube flow and flow between parallel 
plates; however, the results provided in this 
paper are only for parallel plates. 

The pressure gradient term in Eq. (2) is 
assumed to be of oscillatory nature according to: 

t - (S “• 

max 


Assumptions 

Figure 1 shows a sketch of the channel flow 
with the computation domain and boundary condi- 
tions. The following assumptions were made: 


The above equations are subjected to the fol- 
lowing boundary conditions: 


2 


from r = R 

o 

at r = R 0 

at r = 0 

Since the computation domain includes both the 
solid and the fluid, the boundary- condi tion at the 
solid/fluid interface is automatically satisfied. 

It should be noted that the above boundary 
conditions have been carefully selected for the 
fol lowi ng reasons : 

(1) They match the boundary conditions given 
in Kurzweg's [1] analysis. Accordingly, the numer- 
ical results can be compared with his analytical 
solution. 

(2) The above boundary conditions to great 
extent could simulate the central element (from 
both axial and radial directions) of a Stirling 
engine foil regenerator. 



checked against the corresponding analytical veloc- 
ities at several radial locations, figures 2(a), 
(b), and (c) show the numerical velocity versus the 
analytical one at various radial locations in the 
channel (near wall region, y/R 0 = 0.506; half ois- 
tance between the wall and centerline y/R 0 = 0.75; 
and centerline, y/R 0 - 1-0). Results are shown 
for water, w* =* 0.306 (Fig. 2(a)); water, u* = 9.0 
(Fig. 2(b)); and air, w* = 4.48 (Fig. 2(c)). The 
plots show, again the close agreement between the 
numerical and analytical solutions. 

Fluid Flow Resul ts 


Velocity prof i le 

Figures 3(a), (b), and (c) show the dimension- 
less velocity profile u* versus the dimensionless 
distance (y/R 0 ) for water, w* = 0.306 (Fig. 3(a)); 
water, w* = 9 (Fig. 3(b)); and air, oj* = 4.48 
(Fig. 3(c))/ Upon examining the above figures the 
fol lowing were noted: 

(1) For low dimensionless frequency 
(w* = 0.306) the velocity profile has a maximum at 
the channel centerline. 


• Numerical Method 

The above partial differential equations have 
been transformed into finite difference equations 
using different schemes, namely, fully explicit. 
Crank Nicolson and fully implicit. A sensitivity 
analysis was conducted by comparing the results 
from each finite difference scheme with the fluid 
flow analytical solution. It was concluded that 
the Crank Nicolson scheme is the most accurate; 
this is consistent with other research findings 
for unsteady heat equation problems [31. 

RESULTS 

Comparison with Analytical Solution 

One of the main objectives of this work is to 
develop an accurate and efficient numerical scheme 
that can be used with confidence in future engine 
models. Therefore, an extensive grid variation 
test was conducted for both space and time to 
obtain grid independent and accurate numerical 
resul ts . 

Table I shows the normalized tidal displace- 
ment and normalized effective thermal diffusivity 
as obtained for water (Pr = 10.26) with R 0 = 2 R ^ 
and the assumption that the solid has the same 
thermophysical properties as the fluid. In the 
table the present work numerical predictions were 
compared with the analytical solution [1] for w* 
from 0.09 to 9. The numerical results showed the 
same trend where a maximum a* occurs at 
oj* = 0.306 as indicated by the analytical solution 
for Pr = 10.26. Table I shows differences between 
the numerical and analytical results of less than 
2.0 percent. 


(2) As to* increases a boundary layer starts 
to develop near the solid surface. 

(3) Also, as cj* increases the flow direction 
near the wall becomes opposite to that in the core, 
in certain parts of the cycle. 

Friction factor 

From the velocity profiles described above 
one can obtain wall shear stress, friction factor 
and Reynolds number at any pressure gradient phase 
angle. In addition, a root-mean-square can be cal- 
culated for the wall shear stress as well as the 
mean cross-sectional axial velocity. From the lat- 
ter one can calculate friction factor and Reynolds 
number based upon time averaged quantities, simi- 
lar to the work done on circular tubes by Cl>en and 
Griffin [4]. Chen and Griffin had obtained a 
correlation for the normalized friction factor 
(with respect to the corresponding unidirectional 
flow), f/f s , based upon an approximate velocity 
profile given by White C51. Work is underway to 
compare present work results with Chen and 
Griff in' s correlation. 

Heat Transfer Results 

In the channel flow under study, the maximum 
radial temperature difference, between the channel 
wall temperature and the bulk fluid temperature, 
can be obtained by ignoring the radial heat conduc- 
tion. Accordingly, this maximum radial temDera- 
ture difference is equal to the axial temperature 
gradient times the tidal displacement. Now, one 
can obtain a normalized temperature profile utiliz- 
ing the maximum radial temperature difference 
described above. 


Furthermore, the present work velocity pro- 
files for different fluids (water and air) were 


Figures 4(a) and (b) show the normalized tem- 
perature profile versus y/ R 0 for u* = 0.306 
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water), ana a>* = 4.48 (air) at different pressure 
graaient phase angles It should be noted that 
the solid is represented in the region y / R 0 * 0 to 

0.5. r he temperature fluctuation in the fluid is 
much larger, as expected, than it is in the solid. 
Also this ratio between those temperature fluctua- 
tions (in fluid/in solid) becomes greater as w* 
increases . 

Finally, the maximum effective thermal diffu- 
sivity described above in Table I can be better 
explained by plotting both the velocity and temper- 
ature at the same radial location versus 
Figures 5(a) and (b) show the normalized velocity 
and normalized temperature profiles versus <J> at 
y/R 0 * 0.506, 0.75 and 1.0; for . 0.306 
(water), and w* * 4.48 (air) respectively. It 
should be noted that the value of w* chosen for 
each fluid corresponds to the maximum a* that 
can be obtained. " 

It can be seen from the plots that there is 
always a phase shift between the velocity and tem- 
perature profile. It is interesting to notice that 
the phase shift is almost the same for the above 
two cases examined. It is believed that this 
phase shift (approximately 120 ° at y/R 0 * 1.0) 
provides the maximum interaction between axial con- 
vection and radial conduction and consequently the 
maximum net axial heat flux. 

CONCLUDING REMARKS 

An accurate numerical scheme has been applied 
to calculate the fluid flow and heat transfer char- 
acteristics in an oscillating channel flow. The 
flow was assumed laminar, incompressible and has 
constant thermophysical properties. Also a con- 
stant axial temperature gradient was imposed on 
the channel; however the axial heat conduction was 
assumed negligible. 

The numerical scheme has been successful in 
predicting velocity profiles, tidal displacement, 
temperature profiles, and effective thermal diffu- 
sivity for a wide range of «* and different types 
of fluid (water and air). This shows that the 
present work numerical studies confirm earlier 
investigations . 

The interaction between axial convection and 
radial conduction heat transfer, if tuned properly, 
would result in a considerable enhanced axial heat 
transfer. These tuned conditions, for the cases 
examined in the paper, occurred at w* = 0.306 
(water) and 4.48 (air). It is believed that 
the phase shift between the velocity and tempera- 
ture profiles are tied in with this tuning process. 


In a Stirling engine regenerator, these conditions 
(at which maximum axial heat transfer occur) should 
be avoided. 

It should be noted that those tuned u* are 
restricted to the channel configuration and the 
different assumptions explained in the paper. A 
study is underway to examine other operating condi- 
tions (such as type of fluid, wall material, order 
of magnitude of axial pressure and temperature 
gradients, etc.) that are relevant to Stirling 
engines . 

This work is a first step in the planned 
development of a two-dimensional code which can be 
used to model Stirling engine heat exchangers- 
tubular heaters, tubular and foil type regenera- 
tors, and tubular coolers. As the studies 
continue, the simplifying assumptions will be 
changed to include lami nar/turbulent flow, compres- 
sible working fluid, and different outer thermal 
boundary conditions. 

A second study has begun to model the rapid 
expansion and contraction of oscillating flow 
entering and exiting the heat exchangers. This 
work will permit the incorporation of hydrodynamic 
and thermal end effect conditions to the heat 
exchanger model . 
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TABLE I. - COMPARISON BETWEEN PRESENT WORK NUMERICAL RESULTS EOR 
Ax* and a* AND ANALYTICAL SOLUTION [1 ] FOR DIFFERENT u * 
[Water (Pr = 10.26) and R 0 /R j = 2.0.3 


Valensi 

number, 

0)* 

Normalized tidal displacement, 
Ax* 

Normalized effective averaged 
thermal diffusivity, 


Present work 

Analytical [1] 

Present work 

Analytical [1] 

0.09 

0.05941 

0.05995 

0.02616 

0.02568 

.16 

.1055 

.1064 

.03707 

.03645 

.25 

.1643 

.1658 

.04039 

.04029 

.306 

.2019 

.2023 

.04051 

.03979“ 

1 .0 

.6160 

.6181 

.02415 

.02302 

9 

1.607 

1 .61 

.00421 

.00386 




t 



AXIAL 

DIRECTION 


L PARALLEL PLATE 
HEAT EXCHANGER 


(a) SKETCH OF OVERALL SYSTEM. 



'-COMPUTATION 

DOMAIN 


(b) CENTRAL PART SHOWING THE COMPUTATION DOMAIN AND BOUNDARY 
CONDITIONS. 

FIGURE 1. - SKETCH OF A PARALLEL PLATE HEAT EXCHANGER. 
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NORMALIZED VELOCITY U , NUMERICAL SOLUTION 





NORMALIZED VELOCITY U*. ANALYTICAL SOLUTION 
(c) AIR (PT = 0.72), U* = 


FIGURE 2. - COMPARISON BETWEEN NUMERICAL AND ANALYTI- 
CAL VELOCITY PROFILES AT y/R Q = 0.506 (NEAR WALL), 
0.75 AND 1.0 (CENTERLINE) FOR R 0 /Rj = 2.0. (NOTE 
THAT THE DATA POINTS FOR DIFFERENT y/R Q ALMOST 
COINCIDE AND LIE ON A LINE AT 45° WITH THE AXIS). 
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NORMALIZED VELOCITY 


CHANNEL 
CENTERLINE - 


PRESSURE GRADIENT 
PHASE ANGLE , <D 
180 


BSaa*** 








•Wvvvv-yvvw^ 


£ 2222-222 


(a) WATER (Pr = 10.26), R 0 /Rj = 2.0 AND U - 0.306. 









(b) WATER (Pr = 10.26), R Q /R j = 2.0 AND U = 9.0. 



.* V - X **W*X-X*HKX 
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0 .1 .2 .3 .9 .5 .6 .7 .8 

DIMENSIONLESS RADIAL DISTANCE y/R Q 
(C) AIR (Pr = 0.72). R 0 /Rj = 2.0 AND 0* = A. 98. 

FIGURE 3, - NORMALIZED VELOCITY PROFILES VERSUS DIMENSIONLESS RADIAL DISTANCE AT DIFFERENT PRESSURE GRADIENT PHASE 
ANGLES. 







NORMALIZED TEMPERATURE 


CHANNEL 

CENTERLINE 




.8 | — 


SOLID 




N- 


FLUID 


PRESSURE GRADIENT 
PHASE ANGLE, 4> 




(b) AIR CPr = 0.72), R 0 /R j = 2.0 AND U* = A. 48. 

FIGURE 4. - NORMALIZED TEMPERATURE PROFILES VERSUS DIMENSIONLESS RADIAL DISTANCE AT DIFFERENT PRESSURE GRADIENT 
ANGLES. 






NORMALIZED VELOCITY 


VELOCITY 

TEMPERATURE 



(b) AIR (PT = 0.72)/ R 0 /Rj = 2.0 AND U* = A. 48. 

FIGURE 5. - NORMALIZED VELOCITY AND TEMPERATURE PROFILES VERSUS PRESSURE GRADIENT PHASE 
ANGLE AT DIFFERENT RADIAL DISTANCES (y/R Q = 0.506/ 0.75 AND 1.0). 
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NORMALIZED TEMPERATURE 
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